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I. INTRODUCTION 

Heisenberg models are examples of magnetic systems that have true dynamics with the 
real time dynamics governed by coupled equations of motion. According to the classifica- 
tion of the different dynamic universality classes proposed by Hohenberg and Halperin in 
their work on the theory of dynamic critical phenomena [IH, the classical Heisenberg fer- 
romagnet and antiferromagnet are of class J and G, respectively. Although both classes 
have true dynamics, in the former class the order parameter (the uniform magnetization) is 
conserved, whereas in the latter the order parameter (the staggered magnetization) is not a 
conserved quantity. The difference in the dynamic behavior of the Heisenberg ferromagnet 
and antiferromagnet can already be seen from the linear spin-wave theory, which predicts 
different low-temperature dispersion curves for the two models. Dynamic critical behavior is 
describable in terms of a dynamic critical exponent z which depends on conservation laws, 
lattice dimension, and the static critical exponents. The static critical behavior of three- 
dimensional Heisenberg models have been studied using a variety of approaches, including a 
high-resolution Monte Carlo simulation which determined the critical temperature and the 
static critical exponents for simple cubic and body-centered-cubic lattices 0]. In contrast, 
the theory of the dynamic behavior of Heisenberg models is not so well understood. 

A very close realization of an isotropic three-dimensional Heisenberg antiferromagnet is 
RbMnFa. Early experimental studies [0-^ have shown that in RbMnFs the Mn^+ ions, 
with spin S = 5/2, form a simple cubic lattice structure with a nearest- neighbor exchange 
constant of J^^^ = (0.58 ± 0.06) meV and a second-neighbor constant of less than 0.04 meV 
[both defined using our convention for the exchange constant to be shown in Eq. (|T^)]. 



Magnetic ordering with antiferromagnetic alignment of spins occurs below 83K, which we 
denote as the critical temperature T^. In this material, the magnetic anisotropy was found 
to be only about 6 x 10^^ of the exchange field, and no deviation from cubic symmetry 
was seen at Tc 0,01 • Both the static properties and the dynamic response of RbMnFs have 
been examined through neutron scattering experiments. The early work of Tucciarone et 
al ^ found that in the critical region the neutron scattering function has a central peak 
(peak at zero frequency transfer) and a spin-wave peak. Later, the experimental study by 
Cox et al [0] observed a small central peak below Tc as well. The more recent experiments 



by Coldea et al ITDI have also found central peaks for T < T^, in agreement with previous 



experiments. From the theoretical side, renormalization-group (RNG) below T^ |TT]] predicts 
spin-wave peaks, and a central peak in the longitudinal component of the neutron scattering 
function. However, at the critical temperature both renormalization-group [|T2| and mode- 



coupling [r^l theories predict only the spin- wave peak, i.e. the central peak has not been 



predicted by these theories. The central peak is thought to be due to spin diffusion resulting 
from nonlinearities in the dynamical equations 0. In their recent experiment, Coldea et 
al [jlO| obtained the most precise experimental estimate of the dynamic critical exponent, 



z = (1.43 ± 0.04). The theoretical prediction [Q,|T^,|T^] is 2; = 1.5 for class G models in three 
dimensions. 

Large-scale computer simulations using spin-dynamics techniques to study the dynamic 
behavior of Heisenberg ferromagnet and antiferromagnet have been carried out by Chen 



and Landau |14] and Bunker et al |T^, respectively. So far, however, there are no direct 



comparisons of the dispersion curve and the dynamic structure lineshapes obtained from 
simulations with the corresponding experimental results. In the present work we have carried 
out large-scale simulations of the dynamic behavior of the Heisenberg antiferromagnet on 
a simple cubic lattice, and make direct comparisons with experimental data for RbMnFa. 
Sec. II of this paper contains the definition of the model and introduces some simulation 
background. In Sec. Ill we present and discuss our simulation results and compare them 
with experiments. Sec. IV contains some concluding remarks. 



II. MODEL AND METHODS 

A. Model 

The classical Heisenberg antiferromagnetic model is defined by the Hamiltonian 

n=J Yl Sr-Sr', (1) 

<rr'> 

where Sr = (5'r^, Sj.^, Sj.^) is a three-dimensional classical spin of unit length at site r and 
J > is the antiferromagnetic coupling constant between nearest-neighbor pairs of spins. 
We consider L x L x L simple cubic lattices with periodic boundary conditions. The time 
dependence of each spin can be determined from the integration of the equations of motion 



[Eq.(2) in Ref. |]T5|], and the dynamic structure factor S{q,uj) for momentum transfer q and 
frequency transfer u, observable in neutron scattering experiments, is given by 



/+00 f]f 
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where C'^{r — r',t) is the space-displaced, time-displaced spin-spin correlation function de- 
fined, with k = x,y, or z, as 



C\r-r',t) = {Sr'{t)Sr''m - (5/(t))(^r''(0)). (3) 

The displacement r is in units of the lattice unit cell length a. In the case of antiferromagnets, 
the wave-vectors are measured with respect to the (vr, n, n) point which corresponds to the 
Brillouin zone center. 



B. Simulation method 



Using a combination of Monte Carlo and spin-dynamics methods |T^-|1^], we simulated 
the behavior of the simple-cubic classical Heisenberg antiferromagnet with 12 < L < 60 at 
the critical temperature Tc = 1.442929(77) J 0] and below Tc. We have chosen units such 
that the Boltzmann constant ks = ^■ 

Equilibrium configurations were generated using a hybrid Monte Carlo method in which 
a single hybrid Monte Carlo step consisted of two Metropolis steps and eight overrelaxation 
steps |I1|P3[ . Typically 1000 hybrid Monte Carlo steps were used to generate each equihb- 
rium configuration and the coupled equations of motion were then integrated numerically, 
using these states as initial spin configurations. Numerical integrations were performed to 
a maximum time t^ax, using a time step of At. The space-displaced, time-displaced spin- 
spin correlation function C^(r — r', t) was computed for time-displacements ranging from to 
tcutoff- Each of such correlations was calculated from an average over between 40 to 80 differ- 
ent time starting points, evenly spaced by lOAt. Our unit for time is the interaction constant 
J defined in Eq. (|T|). The parameters used in this work were as follows. At T = 0.9Tc we 
used lattices sizes L = 12,24,36,48, and 60, for which the respective time parameters were 
{tmaxJ,tcutoffJ) = (440,400), (440,400), (480,400), (680,600), and (1080, 1000). The respec- 
tive numbers of initial configurations used were 7020, 2850, 1110, 400, and 810. At T^ and 
lattice sizes L = 18, 24, 30, 36, 48 we used tmaxJ = 480 and tcutoff J = 400, and for L = 60 we 
used tmaxJ = 1080 and tcutoff J = 1000. The number of initial configurations used at Tc was 
1000 for L = 18, 30, 36, whereas for L = 24, 48, and 60 the respective numbers were 1125, 715, 
and 510. Other temperatures considered, chosen to coincide with experimental values, were 
T = 0.774Tc, OMGTc and 0.936T,, for which tmaxJ = 480, tcutoff J = 400 and the number of 
initial configurations used was 120. For lattice size L = 24 at T = 0.9Tc the integration was 
carried out with a time step At = O.OIJ"^ using a fourth-order predictor-corrector method. 



as in Refs. p^ , pr5| . For other lattice sizes and temperatures, we used a new algorithm ||T7 
based on fourth-order Suzuki- Trotter decompositions of exponential operators, with a time 
step At = 0.2J^i (except for L = 12 at T = 0.9Tc for which At = O.IJ^^). The latter 
method allowed us to use larger integration time steps and thus to carry out the integration 



to larger tmax, as compared with previous work ||T^,|T5|. In the original study of Chen and 
Landau [14] tmax = 120J~^ and in the work of Bunker et al tmax = 200J~^. In comparison, 
in the present work we have used a larger tmax > 440 J~^, more equilibrium initial spin con- 
figurations, and a larger lattice size, namely L = 60. As it is clear from the above, we have 
concentrated our efforts on two cases: one temperature below Tc, namely T = 0.9Tc, and at 
T 

We also used the technique of calculating partial spin sums "on the fly" |]14|,|15| which 



limits us to data in the (g,0,0), {q,q,0) and {q,q,q) directions with q determined by the 
periodic boundary conditions, 

g=^, n=±l,±2,..,±(L-l),L. (4) 

Since all three Cartesian spatial directions are equivalent by symmetry, the same operation 
carried out for the (g, 0, 0) direction was also carried out for the other two reciprocal lattice 
directions (0, q, 0) and (0, 0, q) and the results were averaged. Similarly, the same operations 
carried out for the {q, q, 0) and (g, q, q) directions were also carried out for the equivalent 
reciprocal lattice directions and the results were averaged for each case. 

For the ferromagnetic case the total magnetization is conserved and the dynamic structure 
factor S'(q, uj) can be separated into a component along the axis of the total magnetization 
(longitudinal component) and a transverse component. However, for the antiferromagnet 
the order parameter is not conserved and separation is not possible. We thus refer to the 
average 

5(q, u) = ^[S^{q, 00) + Sy{q, u) + ^^(q, u)] (5) 

as the dynamic structure factor. 



C. Dynamic finite-size scaling 

Two practical limitations on spin-dynamics techniques imposed by limited computer 
resources are the finite lattice size and the finite evolution time. The finite time cutoff can 
introduce oscillations in S'(q, tu), which can be smoothed out by convoluting the spin-spin 



correlation function with a resolution function in frequency. The finite-size scaling [ll4| , p!8 
can be used to extract the dynamic critical exponent z using 

^^^C^ = G(c.L^gL,5.L^) (6) 

xi(q) 



and 

u^ = L-'Q{qL,6^L'), (7) 

where Sl (q, uj) is the dynamic structure factor convoluted with a Gaussian resolution func- 
tion with characteristic parameter 6^^, COm is a characteristic frequency, defined as 

duj 1 



r^" ^ k, .du i_fc/ N 



and xilq) is the total integrated intensity. 

For the time cutoffs of at least 400 J~^ used in the present simulations, the oscillations in 
the dynamic structure factor due to the finite t cutoff were not very significant. Thus, we first 
estimate the dynamic critical exponent z without using a resolution function, or equivalently, 
we take 5^ = 0. In this case the analysis is simpler and a value for z can be obtained from the 
slope of a graph of log(co'm) vs log(L) (where Um is the characteristic frequency in the absence 
of a resolution function) if the value qL is fixed. It is important to note that lattice sizes 
included in this calculation should be large enough to be in the asymptotic-size regime. The 
approximate lattice size of the onset of the asymptotic regime can be estimated by looking 
at the behavior of ujL^ for different lattice sizes using trial values of z. 

The effects of the small oscillations in the dynamic structure factor on the dynamic 
critical exponent z can be evaluated by repeating the analysis using a resolution function. 
For this purpose, we chose 

"60 



S... = 0.02 



(9) 



L 
so that the function Q{qL, 6^L^) in Eq. (^ is a constant if qL is fixed, yielding 

U^m ~ L-\ (10) 

Because (5^ depends on z, this exponent had to be determined iteratively. We used n = 1,2 
and several initial values for z in the iterations, in order to check how stable the converged 
value of z is. 

III. RESULTS 

A. Numerical data for 5(q, lo) 

For T < Tc our results for the dynamic structure factor, as defined in Eq. (^, show a spin- 
wave and a central peak. In Fig. |l] we show lineshapes for lattice size L = 60 at T = 0.9Tc 



and wave-vectors q = vr/lO, 7r/6 and vr/S in the [100] direction. We see that as q increases, 
the central peak broadens and its relative amplitude increases. Fig. ^ shows lineshapes for 
L = 60 at Tc and wave- vectors q = vr/lO and 7r/6 in the [100] direction. It is clear from these 
lineshapes that the oscillations due to the finite tcutoff are indeed negligible; therefore, in our 
analysis of the lineshapes we have not convoluted our results with a resolution function. (As 
explained in the following section, we later convoluted our results with a Gaussian resolution 
function in order to directly compare our lineshapes with the experiments. The reason for 
this is that there is an intrinsic finite resolution in the experimental data due to the finite 
divergence of neutron beams.) The width of any structure in the lineshapes discussed here 
is much larger than our resolution in frequency Aa; = 1.2tt /tcutoff ■ 

Below Tc, previous theoretical |jll[ and experimental ||10[ studies motivated us to extract 



the position and the half-width of the spin-wave and central peaks by fitting the lineshape 
to a Lorentzian form 

A^l BTl BTl 

where the first term corresponds to the central peak and the last two terms are contributions 
from the spin-wave creation and annihilation peaks at cu = ±ujs. For T = 0.9Tc we find that 
Lorentzian lineshapes fit our results reasonably well for small values of q, except for the 
smallest value, namely q = 27t/L, in the [100] direction. The reason is that for each lattice 
size L, the dynamic structure factor for the smallest value of q corresponds to correlations 
between spins displaced in space by a distance L/2, and the effect of the finite lattice size 
is particularly prominent in these cases, causing the lineshapes to depart significantly from 
a Lorentzian form. For large values of q (approximately q > 27r(L/4 — 2)/L) the Lorentzian 



form given in Eq. {\iT\ ) does not fit the data, especially at large frequency transfer. In 
general, the fitted parameters varied when different frequency ranges were used in the fit. 
Although this variation was small, it was often larger than the statistical error in the fitted 
parameters obtained from the fit using a single frequency range. Therefore, for T = O.QT^ 
we estimated the error in the fitted parameters by fitting the lineshapes using three different 
ranges of frequency. The values of the parameters were then averaged and the error bars 
estimated from the fluctuations. At Tc, renormalizat ion-group theory (RNG) |T^ predicts a 
non-Lorentzian functional form for the spin-wave lineshape, which has been used along with 
a Lorentzian central peak to analyze experimental data [|10|- Since it is more complicated 



to perform fits to this RNG functional form and since the spin-wave peaks obtained from 
the simulations are more pronounced than in the experiment, and thus less dependent on 
the fitted functional, we have fitted the lineshapes at Tc to Lorentzians, as given in Eq. 
(pA|). Although obtaining a good fit to our data at Tc was more difficult than below Tc, the 
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resulting fits at Tc are still reasonable. Unlike for T = 0.9Tc, at Tc the lineshape parameters 
used in the analysis below are the values obtained from the fit to only one frequency range, 
which was the one that gave the best fit. One should then expect that the actual error in 
the fitted parameters at Tc is larger (by up to a factor of 5) than the error bars shown in the 
figures below. 

In addition to the spin-wave and the central peaks, we observed some other peaks on 
the large frequency tail of the spin-wave peaks. Although these large frequency peaks had 
very small amplitudes, they could be discerned from the background fluctuations. Using 
the spin-wave frequencies in the [100], [110] and [111] directions we could check that the 
position of these extra peaks corresponded to frequencies of two spin-wave addition peaks. 
These extra structures in the lineshapes were particularly visible for the smallest values of 
wave- vectors. 

Fig. ^ shows how the dispersion curve varies as the temperature increases from T = 
0.774Tc to Tc. The dispersion curves illustrated here are for the [100] direction and are 
plotted up to g = 7r/2, which corresponds to one half of the Brillouin zone. As mentioned 
before, for larger values of q the Lorentzian in Eq. ( pT]) did not yield good fits to the 
lineshapes; however, the spin-wave positions could still be directly read off the graphs, but 
with larger error bars. For our present purpose of observing the approach to the critical 
region as the temperature is raised from below Tc, it suffices to consider the dispersion curve 
for wave-vectors up to g = 7r/2. Well below Tc, the dispersion relation is linear for small q; 
as T — > Tc, it changes gradually from linear to a power-law behavior of the form 

tOs = A,q\ (12) 

For T = 0.9Tc and L = 60 a fit including the smallest five q values of the dispersion curve 
to Eq. (|1^) yielded x = 1.017 ± 0.003. As we probed further away from the Brillouin zone 
center by including larger values of q in the fit, the exponent decreased slightly. In order to 
check how sensitive the fitted exponent is to the particular form of the fitted function, we 
have performed new fits to a function which includes a quadratic term, i.e.. 

Us = A.g^ + S,g2. (13) 

Fitting the smallest five q values of the dispersion curve to a function of the form given 
by Eq. (|13D yielded an exponent x = 1.020 ± 0.003, which is in good agreement with the 
value obtained from the previous fit. When larger values of q in the dispersion curve were 
included in the fits, Eq. (|l^) tended to yield smaller x^'s per degree of freedom than Eq. 



(|T^). The dispersion curve for T = Tc and L = 60 fitted to Eq. ([1^ ) yielded an exponent 



of a: = 1.38 ± 0.01 when the smallest 12 values of q were included in the fit. As the larger 
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q were excluded from the fit, the exponent increased slightly, tending towards x ~ 1.40. 
When only the smallest few values of q were included in the fit, the exponent decreased 
again, reflecting the fact that as we probed correlations between spins separated by larger 
distances (or equivalently, smaller q) the finite size of the lattice (and thus of the correlation 
length) is revealed, showing that the system is not at criticality. Hence, the exponent x 
decreases towards unity. On the other hand, large values of q correspond to short distance 
(in the direct lattice space) spin-spin correlations, and the correlation length is much larger 
than the distance probed. One would thus expect the critical behavior of the system to 
be manifest, and this is indeed consistent with what we obtained. Our results at T^. are in 



agreement with the recent experiment [T^ which obtained an exponent x = 1.43 ±0.04 when 
the dispersion curve at T^ was fitted to a power-law relation of the form given in Eq. ([T^) . 
The solid lines in Fig. ^ are fits to Eq. (|13D; in general, these fits gave lower values of x^ 
per degree of freedom than fits to Eq. (|12D. 

In the critical region, dynamic scaling theory predicts that the half-width of spin-wave 
peaks behaves as a power-law, r2 ~ q^'^ H^, whereas for the hydrodynamic regime the 
prediction from hydrodynamic theory is r2 ~ q^ [|3l- The half- width of the spin- wave 
peaks at T = O.QT^ and L = 60 from our simulations is shown in Fig. ^. We observed 
a crossover from r2 = (0.401 ± 0.004)g^'^^^°°^ for larger values of wave-vector to r2 = 
(0.48 ± 0.02)q^'^^^^'^^ for small values of q. The behavior for relatively large wave-vectors is 
in agreement with dynamic scaling theory and with the recent experiment |T^ . The exponent 
we obtained by fitting only small values of q is close to the hydrodynamic prediction. Thus, 
for the spin-wave half-width we have observed a crossover of exponents associated with two 
different regimes, namely, the critical and the hydrodynamic regions. This crossover is similar 
to the one observed in the dispersion curve at Tc, discussed above. For T = Tc and L = 60 
the spin-wave half-width also had a power-law behavior which varied from approximately 
g^'^ when the 12 smallest values of q were included to approximately q^'^ when only the 
smallest five wave- vectors were included in the fit. In their recent experiment, Coldea et al 
T0| obtained r2 = /}gi-4i±o-05 f^j, ^j-^g temperature range 0.77Tc < T < Tc, and the coefficient 



D increased with increasing temperature. 

As in the experiments, the dynamic structure factors from our simulations had central 
peaks (zero frequency transfer peaks) for T < T^. In contrast, renormalization-group theory 
predicts a central peak in the longitudinal component of the dynamic structure factor only 
below Tc []TT|, and none of the theories predict a central peak at T^ [|T^,|T^. For T = O.QT^ 
and L = 60 fitting the central peak half-width to the form Fi ~ q^ yielded very large x^ per 
degree of freedom. A much improved fit was obtained by using the function Fi = Ai + Biq^^ , 
which allows for a non-zero central peak width when q vanishes. In these fits the data for the 



smallest value of q observable in L = 60, i.e. n = 1, were not included, because of the large 
finite-size effects in them. The fit including data for q corresponding to n = 2 until n = 7 
yielded Ai ~ 0.013 ± 0.001, Bi ~ 0.120 ± 0.005 and Ci ~ 2.4 ± 0.2. As we systematically 
included larger values of q in the fits, these parameters decreased slightly. At Tc we have 
also fitted the central peaks to Lorentzians, according to Eq. (0); however, these tended to 
yield curves with smaller amplitudes than the data, as can be seen in Fig. ^j. Since there is 
no theoretical prediction for the central peak, we have also tried to fit them with a Gaussian 
form. These latter fits were, nevertheless, much worse than the fits with Lorentzians. 

The lattice sizes that we used, namely L = 12, 24, 36, 48, and 60, are all multiples 
of 12; thus, there are certain wave- vectors which are common to all lattice sizes. This is 
an advantage in the study of effects due to finite lattices, because it allows us to compare 
lineshapes and spin-wave dispersion relations for different lattice sizes at a fixed value of 
wave-vector. At T = 0.9Tc we did not see a significant finite-size effect for L > 24; however, 
when we superimposed lineshapes at Tc for a fixed value of q, and different values of L, 
finite-size effects were noticeable for L = 24. For the larger values of L that we used, the 
lineshapes were the same within the error bars. 

The dynamic critical exponent z was extracted from the finite-size scaling of u!m, as 
described in a previous section. We started the analysis using no resolution function, or 
equivalently 6uj = 0, and n = 1,2. As in previous work [0, we estimated the lattice L = 30 
to be approximately the onset of the asymptotic-size regime. From the slope of a \og{uJm) 
vs log(L) graph fitted with four data points, namely lattice sizes L = 30, 36, 48, and 60, we 
obtained z = 1.45 ± 0.01 for n = 1 and z = 1.42 ± 0.01 for n = 2. In order to check the 
effects of the very small oscillations in S{q, u) due to the finite cutoff time, we proceeded 
to estimate the value of z using a resolution function with 6^ given by Eq. @. For the 
iterations, we used several initial values of z^'^\ ranging from z^^'^ = 1.31 to 1.59, and in 
all cases the exponent z converged rapidly to the final value with at most three iterations 
being necessary. For all the initial values of z^'^^ that we used, we obtained an exponent 
z = 1.43 ± 0.01 for n = 1 and z = 1.42 ± 0.01 for n = 2. In general, the x^ P^r degree of 
freedom in the fits for n = 2 were slightly lower than for n = 1, although it was reasonable 
in all cases. Our final estimate for the dynamic critical exponent is z = 1.43 ± 0.03, where 
the error bar reflects the fluctuations in the different estimates of z. A comparison of the 
characteristic frequency Um as a function of the lattice size L for the analysis with and 
without a resolution function is shown in Fig. |^. For the former case, the graph shown 
corresponds to the converged values of z for both n = 1 and 2. 



B. Comparison ^vith experiment 

In this section we compare our results with the recent neutron scattering experiment by 
Coldea et al |[1^ . Before proceeding with the direct comparison, it is necessary to clarify the 
units and possible normalization factors between simulation and experiment. 

The neutron scattering experiment was done on RbMnFs, which can be described by a 
quantum Heisenberg Hamiltonian of the form 

n = r^'p Y. Sr"^ • Sr''^, (14) 

<rr'> 

where Sr are spin operators with magnitude |Sr P = S{S + 1) and the interaction 
strength between pairs of nearest-neighbors was determined experimentally to be J^^^ = 
(0.58 ± 0.06) meV . In contrast, our simulations were performed on a classical Heisenberg 
Hamiltonian given in Eq. {^. However, quantum Heisenberg systems with large spin values 
{S > 2) have been shown to behave as classical Heisenberg systems, where the spins are taken 



to be vectors of magnitude \/S{S + 1) with the same interaction strength between pairs of 



nearest neighbors as in the quantum system |21|. In our simulations the spins were taken 



to be vectors of unit length. Hence, to preserve the Hamiltonian the interaction strength J 
from our simulation has to be normalized according to 

J = J--P5(5 + l). (15) 

Although this choice of normalization for spin vectors and the interaction strength leaves 
the Hamiltonian unchanged, it does modify the equations of motions. The dynamics of the 
classical system so defined is the same as the quantum system defined by the Hamiltonian 
in Eq. (P^ if one rescales the time, or equivalently, the frequency transfer. We obtain 



where u^^^ is the frequency transfer in the quantum system, measured experimentally, and 
w/J is the frequency transfer in units of J from our simulations. Parenthetically, we note 
that the critical temperature of the classical Hamiltonian in Eq. (|I]) has been determined 
from simulations to be Tc = 1.442929(77) J 0. Using the normalization for the interaction 
strength J given in Eq. ([T5|) and the experimental value J^^^ = (6.8 ± 0.6)K [^ we get 
Tc = (85.9 ± 7.6)K, where the 9 percent error comes from the uncertainty in J'^^^. The 
experimental value of the critical temperature is around 83K which is well within the error 
bars. 

Due to detailed balance, neutron scattering experiments measure the dynamic structure 
factor multiplied by a temperature and frequency dependent population factor [jlO|,^^ . 
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This factor does not appear in the simulations of the classical system for which the dynamic 
structure factor is computed directly. For the comparison, we removed the population factor 
from the experimental data. 

The finite divergence of the neutron beam gives rise to a resolution function which is 
usually approximated by a Gaussian in the 4-dimensional energy and wave- vector space. 
In the experiment [|l^, the measured resolution width along the energy axis was 0.25 meV 



(full-width at half-maximum) for incoherent elastic scattering. In order to directly compare 
our results with the experiment, we convoluted the lineshapes from our simulation with a 
Gaussian resolution function in energy with the experimental value of full-width at half- 
maximum, normalized according to Eq. ([I6|). The standard deviation 5^ thus obtained 



for the Gaussian resolution function was 0.0619 in units of J . As a test of the effects of 
the resolution function in the wave-vector space, we have also convoluted our lineshapes 
with a 3-dimensional Gaussian function where, for simplicity, we have taken the resolution 
width in the three wave-vector components to be the same, and equal to the average of 
the experimental resolution in the longitudinal and transverse directions. The effect of this 
convolution was found to be negligible; thus the lineshapes used in the comparisons shown 
below do not include the resolution in wave-vector. 

The experiment |l^ performed constant wave-vector scans with both positive and nega- 



tive energy transfer. The wave-vector transfer Q was measured along the [1,1;1] direction, 
around the antiferromagnetic zone center which in our notation is the (vr, tt, tt) point. Note 
that Ref. |]10| defines the wave-vector transfer Q in units such that the antiferromagnetic 
zone center is (0.5,0.5,0.5); hence, to express Q of Ref. flOl in units of A~^ one has to 



multiply it by 27r/a, where a is the cubic lattice parameter expressed in A. However, in 
the simulation direct lattice positions are defined in units of the lattice constant a; thus we 
obtain wave-vectors multiplied by the constant a. Let us emphasize that one has to divide 
the wave-vector Q [and also g, see Eq. (H)] defined in this paper by 2ti in order to express it 
in the units used in Ref. |]10|- In the experiment, measurements were taken for wave- vectors 
Q = (yr+g, TT+g, 7r+g), with q = 27r(0.02), 27r(0.04),..., 27r(0.12), but unfortunately these val- 
ues of q are not all observed in our simulations with the particular lattice sizes that we used. 
For instance, with a lattice size L = 60 we observe wave-vectors with q = 27r(0. 01666...), 
27r(0. 03333. ..),... and so on, according to Eq. (0). Thus, in order to directly compare the 
lineshapes from the simulation with the experimental ones, it was necessary to interpolate 
our results to obtain the same q values of the experiment. This was done by first fitting our 



lineshapes with a Lorentzian form, as given in Eq. ([TlD . Since the parameters B, T2 and Ug 
obtained from these fits behave as power-laws of g, we linearly interpolated the logarithm 
of these parameters as a function of the logarithm of g, to obtain new parameters for the 
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lineshapes corresponding to those values of q observed in the experiment. We estimated the 
uncertainties from this procedure to be less than five percent for the parameter B, less than 
three percent for the spin-wave half-width r2 and the spin-wave position ujs at Tc-, and less 
than one percent for the spin-wave position ujs aX T = 0.9Tc. Below T^., the parameters A 
and Fi associated with the central peak were linearly interpolated, yielding new parameters 
with uncertainties of approximately five percent. At Tc, the parameter A was interpolated 
in the log-log plane (as for B, V2 and ojs discussed above), whereas Fi was simply linearly 
interpolated. The uncertainties in A and Fi at T^. were estimated to be less than ten per- 
cent. For L = 60, there is one value of q, namely q = 27r(0.10), which is accessible to both 
simulation and experiment. This was the only case for which we did not have to interpolate 
in q. 

Below Tc, the simulations are mainly for T = 0.9Tc which unfortunately does not coincide 
with any temperature used in the experiment; however, it is very close to T = 0.894Tc which 
is one of the temperatures for which experimental results are available. To correct for the 
slight difference, we made a linear interpolation in temperature, using our results for L = 24 
at T = 0.846Tc and at T = 0.9Tc. We first fitted the lineshapes at these two temperatures 
to a Lorentzian of the form given by Eq. ([lID , then we linearly interpolated the position and 
the amplitude of the spin- wave peak at these temperatures, to obtain the spin- wave position 
and amplitude corresponding to T = 0.894Tc. For small values of q we found that the 
frequency of the spin-wave peak at T = 0.894Tc was approximately 1.5 percent larger than 
at T = 0.9Tc and this difference decreased for larger values of q. The spin-wave amplitude 
at T = 0.894Tc was found to be approximately 5 percent larger than at T = 0.9Tc for small 
values of q. As for the spin-wave position, the difference in the amplitudes decreased for 
larger values of wave-vector. 

The intensity of the lineshapes in the neutron scattering experiment was measured in 
counts per 15 seconds. For both temperatures T = 0.894Tc and T = Tc the measurements 
for the several wave-vectors were done with the same experimental set-up and conditions. 
Therefore, the relative intensities of the lineshapes for the different wave-vectors is fixed, 
and equal for both temperatures. The intensity of the lineshapes obtained in the simulation 
had to be normalized to the experimental value; however, because the relative intensities 
for different wave-vectors is fixed, we have only one independent normalization factor for 
all the wave-vectors at both temperatures. The normalization of the intensity was chosen 
so that the spin-wave peak for T = 0.894Tc and q = 27r(0.08) from the experiment and the 
simulation matched. This same factor was used to normalize the intensities of the lineshapes 
corresponding to the remaining values of wave-vectors at T = 0.894Tc, and for all cases at 
Tc. The factor used was 70 counts/15 sees., which we multiplied to the simulated lineshapes 
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for all values of g at T = 0.894Tc and at Tc. 

The final lineshapes for T = 0.894Tc, L = 60, and several wave-vectors are shown in Fig. 
P together with experimental lineshapes for each case. Figs. |^(a) and 0(b) show respectively 
the comparisons of the dispersion curve and the spin-wave half-width from the simulation 
and the experiment at T = 0.894Tc. The good agreement between our results and experiment 
can be seen from either the direct comparison of the lineshapes, or the comparisons of the 
dispersion curve and the spin-wave half-width. There is an agreement between the lineshape 
intensities from simulation and experiment over two orders of magnitude, from q = 27r(0.02) 
to g = 27r(0.10). Fig. ^ shows the comparison of lineshapes from the simulation and the 
experiment for T = T^, L = 60, and several values of q. The dispersion curve obtained from 
the simulations at T = T^, shown in Fig. P, is systematically larger than the experimental 
values. We would like to emphasize that the error bars shown for the dispersion curve 
obtained from our simulations at Tc reflect only the statistical errors of a best fit of the 
lineshapes with Eq. (11). For each wave- vector, this fit was done with only one range of 



frequency; hence errors associated with the choice of frequency range and the quality of 
the fit were not taken into account. It is reasonable to expect that such sources of error 
would increase the error bars by a factor of 5. From the direct comparison of the simulated 
and experimental lineshapes at Tc it is difficult to determine the difference in the spin-wave 
frequencies, because the spin-wave peaks from the experiment are not very pronounced, 
and their positions have to be extracted from the fits of the lineshapes. As we mentioned 
before, the experimental data at Tc were fitted to a functional form predicted by RNG 
theory plus a Lorentzian central peak. As an illustration, one such fit is shown in Fig. H(c) 
for q = 27r(0.08), along with the RNG component of the fit and the prediction by mode- 
coupling theory. Finally, even though at Tc the lineshape intensities from the simulations for 
small frequency transfer tended to be lower as compared to the experiment, the agreement 
is still reasonably good, considering the variation of the intensities over almost two orders of 
magnitude from q = 27r(0.02) to g = 27r(0.12). 



IV. CONCLUSION 

We have studied the dynamic critical properties of the classical Heisenberg antiferromag- 
net in a simple cubic lattice, using large-scale computer simulations. A new time integration 
technique implemented by Krech et al allowed us to use a larger time integration step and 
we were thus able to extend the maximum integration time to much larger values than in 
previous work. 
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Below Tc, the dispersion curves were approximately linear for wave-vectors well within 
the first Brillouin zone. Increasing the temperature towards the critical temperature the 
dispersion curve became a power-law, reflecting the crossover from hydrodynamic to critical 
behavior. The power-law behavior of the spin-wave half-width at T = 0.9Tc also showed 
a crossover from critical behavior at large values of q to hydrodynamic behavior at small 
values of q. The dynamic critical exponent was estimated to be 2 = (1.43 ± 0.03) which is in 
agreement with the experimental value of Coldea et al, and slightly lower than the dynamic 
scahng prediction. 

We made direct, quantitative comparison of both the dispersion curve and the lineshapes 
obtained from our simulations with the recent experimental results by Coldea et al. At 
T = 0.894Tc the agreement was very good. The major difference was at T^ where spin-wave 
peaks from our simulations tended to be at slightly larger frequencies than the experimental 
results. Both at T = 0.894Tc and at Tc the lineshape intensities varied over almost two 
orders of magnitude from q = 27r(0.02) to g = 27r(0.10) and there was good agreement 
between the intensities from simulation and experiment over the whole range. Thus, the 
simple isotropic nearest-neighbor classical Heisenberg model is very good for describing the 
dynamic behavior of this real magnetic system, except for small differences in spin-wave 
frequencies at the critical temperature. 
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FIGURE CAPTIONS 

Fig. 1: Dynamic structure factor ^(q, u) from our simulations for L = 60 at T = 0.9Tc 
and wave-vectors (a) q = tt/IO, (b) tt/6 and (c) n/3 in the [100] direction. The symbols 
represent spin dynamics data and the solid line is a fit with the Lorentzian function 
given in Eq. ([TT|). For clarity, error bars are only shown for a few typical points, i.e. 
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error bars for the data in the neighborhood of each of these points are similar. At high 
frequencies error bars are of the size of the fluctuations in these data. 

Fig. 2: Dynamic structure factor ^(q, uj) from our simulations for L = 60 at T = Tc and 
wave-vectors (a) q = vr/lO and (b) tt/6 in the [100] direction. The symbols represent 
spin dynamics data and the solid line is a fit with the Lorentzian function given in Eq. 
11). As in Fig. |l], error bars are only shown for a few typical points. 



Fig. 3: Spin-wave dispersion relations for T < T^ in the [100] direction. The symbols 
represent spin-wave positions extracted from Lorentzian fits to the lineshapes from 
the simulations, and the solid curves are fits of the dispersion relations at different 



temperatures to Eq. (plsl). 



Fig. 4: Log-log graph of the half- width of the spin-wave peak extracted from Lorentzian 
fits to the lineshapes obtained from simulations for L = 60 and T = 0.9Tc in the [100] 
direction as a function of q. 

Fig. 5: Finite-size scaling plot for ujm (with qL =const, 6i^L^ =const) for the analysis 
with and without a resolution function. For the former case, the data used correspond 
to the converged values of z, for n = 1,2. The error bars were smaller than the symbol 
sizes. 

Fig. 6: Comparison of lineshapes obtained from fits to simulation data for L = 60 
(solid line) and experiment (open circles) at T = 0.894Tc in the [111] direction: (a) 
q = 27r(0.04), (b) q = 27r(0.06), (c) q = 27r(0.08), and (d) q = 27r(0.10). The horizontal 
line segment in each graph represents the resolution in energy (full-width at half- 
maximum) . 

Fig. 7: Comparison of the (a) dispersion curve and the (b) spin-wave half-width, 
obtained from simulations for L = 60 (open circle) and the experiment (open triangle) 
at T = 0.894Tc, in the [111] direction. The simulation data shown here correspond to 
values of q accessible with L = 60, without interpolation to match the q values from 
the experiment. 

Fig. 8: Comparison of lineshapes obtained from fits to simulation data for L = 60 (solid 
line) and experiment (open circles) at T = T^ in the [111] direction: (a) q = 27r(0.04), 
(b) q = 27r(0.06), (c) q = 27r(0.08), and (d) q = 27r(0.10). The dot-dashed line in (c) 
is a fit of the experimental data to the functional form predicted by the RNG theory 
plus a Lorentzian central peak, and the RNG component of the fit is shown by the 

15 



long-dashed line. The prediction of Mode Coupling (MC) theory for q = 27r(0.08) is 
shown by the dotted line in (c). The horizontal line segment in each graph represents 
the resolution in energy (full-width at half-maximum). 

Fig. 9: Comparison of the dispersion curve obtained from our simulation for L = 60 
(open circle) and the experiment (open triangle) at T = T^, in the [111] direction. In 
the notation here, the first Brillouin zone edge is at \{q,q,q)\ — 2.72. The simulation 
data shown here correspond to values of q accessible with L = 60, without interpolation 
to match the q values from the experiment. 
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